% Load the data file
load('D:\Project\Data\carbon_with2point5mmbone\carbon_with2point5mmbone\carbon_with1point5mmbone.mat');

%% Basic Parameter Setup
% These parameters are defined by the experimental setup and are consistent
% for both PA and US imaging.

% General acquisition parameters
params.fs = Receive(1).decimSampleRate * 1e6; % [Hz] Sampling frequency
params.fc = Trans.frequency * 1e6;          % [Hz] Transducer center frequency
params.c = 1500;                            % [m/s] Speed of sound in the medium

% Transducer geometry
params.Nelements = 64;
params.pitch = 0.300e-3;    % [m] Element pitch
params.width = 0.280e-3;    % [m] Element width
params.radius = inf;        % [m] Linear array

% Imaging parameters
params.bandwidth = 90;      % [%] Relative bandwidth
params.fnumber = [];        % Use automatic f-number calculation in DAS

%% Data Separation (PA and US)
% This section separates the interleaved PA and US data from the raw RcvData.

% Define the number of acquisitions per frame for PA and US
ne = 5; % Number of Photoacoustic (PA) acquisitions per frame
na = 5; % Number of Ultrasound (US) acquisitions per frame

% --- Photoacoustic (PA) Data Extraction ---
RcvData_PA = {RcvData{1}(1:Receive(ne).endSample,:,:)};

% Create indices to select only the PA acquisitions
Indices1_PA = repelem((0:(size(RcvData_PA{1},3)-1)).*(ne+na), ne);
Indices2_PA = repmat(1:ne, 1, size(RcvData_PA{1},3));
PAReceiveIndices = Indices1_PA + Indices2_PA;
Receive_PA = Receive(PAReceiveIndices);

% Reorder frames to match acquisition sequence
frame_order = [Resource.RcvBuffer.lastFrame+1:Resource.RcvBuffer.numFrames, 1:Resource.RcvBuffer.lastFrame];

% Read and organize PA data
data_PA = zeros(numel(Receive_PA(1).startSample:Receive_PA(1).endSample), Trans.numelements, ne, numel(frame_order));
n_pa = 1;
frame_idx_pa = 0;
for n_frame_pa = frame_order
    frame_idx_pa = frame_idx_pa + 1;
    for n_tx_pa = 1:ne
        data_PA(:,:,n_tx_pa,frame_idx_pa) = RcvData_PA{1}(Receive_PA(n_pa).startSample:Receive_PA(n_pa).endSample, Trans.Connector, n_frame_pa);
        n_pa = n_pa + 1;
    end
end

% --- Ultrasound (US) Data Extraction ---
RcvData_US = {RcvData{1}(Receive(ne+1).startSample:Receive(ne+na).endSample,:,:)};

% Create indices to select only the US acquisitions
Indices1_US = repelem((0:(size(RcvData_US{1},3)-1)).*(ne+na), na);
Indices2_US = repmat((ne+1):(ne+na), 1, size(RcvData_US{1},3));
USReceiveIndices = Indices1_US + Indices2_US;
Receive_US = Receive(USReceiveIndices);

% Adjust sample indices for the extracted US data block
sampleDisplacement = Receive_US(1).startSample - 1;
for i = 1:numel(Receive_US)
    Receive_US(i).startSample = Receive_US(i).startSample - sampleDisplacement;
    Receive_US(i).endSample = Receive_US(i).endSample - sampleDisplacement;
end

% Read and organize US data
data_US = zeros(numel(Receive_US(1).startSample:Receive_US(1).endSample), Trans.numelements, na, numel(frame_order));
n_us = 1;
frame_idx_us = 0;
for n_frame_us = frame_order
    frame_idx_us = frame_idx_us + 1;
    for n_tx_us = 1:na
        data_US(:,:,n_tx_us,frame_idx_us) = RcvData_US{1}(Receive_US(n_us).startSample:Receive_US(n_us).endSample, Trans.Connector, n_frame_us);
        n_us = n_us + 1;
    end
end

%% Image Reconstruction
% This section performs beamforming and creates the images.

% --- IQ Demodulation ---
% Convert RF data to IQ (in-phase and quadrature) data.
% The initial samples are removed to exclude laser trigger delays or transmit artifacts.
IQ_PA = rf2iq(data_PA(9:end,:,1,1), params.fs, params.fc);
IQ_US = rf2iq(data_US(8:end,:,1,1), params.fs, params.fc);

% --- Imaging Grid Setup ---
xlen = params.pitch * (params.Nelements - 1) + params.width;
x_imaging_axis = (-xlen/2):1e-4:(xlen/2);
z_imaging_axis = 0:0.1e-3:40e-3;
[xi, zi] = meshgrid(x_imaging_axis, z_imaging_axis);

% --- PA Reconstruction ---
params_PA = params;
params_PA.passive = true; % PA is a passive modality (receive-only)
bfsig_PA = das(IQ_PA, xi, zi, zeros(1, params.Nelements), params_PA);

% --- US Reconstruction ---
params_US = params;
bfsig_US = das(IQ_US, xi, zi, zeros(1, params.Nelements), params_US);


%% Visualization
% This section displays the reconstructed PA and US images.

% --- Display PA Image ---
log_compressed_parameter_pa = 16;
figure;
IPA = bmode(bfsig_PA, log_compressed_parameter_pa);
imagesc(x_imaging_axis * 1000, z_imaging_axis * 1000, IPA);
colormap(hot);
title('PA Image (DAS)', 'FontSize', 16);
xlabel('Lateral Distance (mm)');
ylabel('Depth (mm)');
colorbar;
axis image;

% --- Display US Image ---
log_compressed_parameter_us = 40;
figure;
IUS = bmode(bfsig_US, log_compressed_parameter_us);
imagesc(x_imaging_axis * 1000, z_imaging_axis * 1000, IUS);
colormap(gray);
title('US Image (DAS)', 'FontSize', 16);
xlabel('Lateral Distance (mm)');
ylabel('Depth (mm)');
colorbar;
axis image;

% --- Display Zoomed Images ---
figure;
imagesc(x_imaging_axis * 1000, z_imaging_axis(201:end) * 1000, IPA(201:end,:));
colormap(hot);
title('Zoomed PA Image', 'FontSize', 16);
xlabel('Lateral Distance (mm)');
ylabel('Depth (mm)');
axis image;

figure;
imagesc(x_imaging_axis * 1000, z_imaging_axis(201:end) * 1000, IUS(201:end,:));
colormap(gray);
title('Zoomed US Image', 'FontSize', 16);
xlabel('Lateral Distance (mm)');
ylabel('Depth (mm)');
axis image;


%% Profile Analysis
% Compare the signal profile from the PA image with a reference.
load('D:\Project\Data\carbon_with2point5mmbone\carbon_with2point5mmbone\cc.mat'); % Load reference profile data

figure;
plot(cc, 'linewidth', 2);
hold on;
plot(IPA(343,:), 'linewidth', 2);
title('Effect of Bone on Target Imaging');
legend({'Target without bone', 'Target with bone plate'});
grid on;
xlabel('Lateral Position');
ylabel('Amplitude');
set(gca, 'FontSize', 12);